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ABSTRACT 

main  result  of  this  paper  is  a  stability  theorem  for  a  certain  class 
of  difference  algorithms  designed  to  give  approximate  solutions  of  a  model 
inverse  scattering  problem  in  one  dimension.  This  stability  result  guarantees 
the  convergence  of  the  approximate  solutions  to  the  exact  solution  of  the 
problem  as  the  grid  of  the  difference  scheme  is  refined.  ■Wa  prasoilt  the 
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results  of  numerical  experiments^Dased  on  one  of  these  schemes,  in  which  second- 
order  convergence  is  observed.  Furthermore  the  cost  (that  is,  the  dependence 
on  N  of  the  number  of  arithmetic  operations  required  to  compute  the  solution 
at  N  grid  points)  of  the  algorithms  discussed  below  is  essentially  optimal. 
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SIGNIFICANCE  AND  EXPLANATION 


Many  problems  of  applied  mathematics  involve  the  inference  of  mechanical 
properties  of  a  medium,  parts  of  which  are  not  accessible  to  direct  observation, 
from  measurements  of  scattering  of  small-amplitude  waves.  Such  problems 
arise  in  exploration  geophysics,  physical  chemistry,  and  ultrasound  tomography, 
among  other  areas. 

Many  of  these  problems  are  equivalent  in  principle  to  boundary  value 
problems  for  certain  partial  differential  equations.  Any  method  of  approximate 
solution  for  such  boundary  value  problems  must  have  the  attribute  of  numerical 
stability,  in  order  to  be  useful:  that  is,  the  round-off  errors  present  in 
all  computation  must  not  lead  to  uncontrolled  growth  in  errors  of  the  computed 
quantities. 

In  this  paper,  a  class  of  difference  algorithms  for  a  simple  model 
inverse  scattering  problem  is  shown  to  be  stable.  This  result  implies  that 
these  algorithms  will  actually  compute  approximate  solutions  to  the  inverse 
scattering  problem.  In  fact,  explicit  error  bounds  can,  in  principle,  be 
extracted  from  the  results  presented  here. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC ,  and  not  with  the  author  of  this  report. 
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NUMER1CAL  STABILITY  IN  AN  INVERSE 
SCATTERING  PROBLEM 


W.  W.  Symes 


51.  Introduction 

The  main  result  of  this  paper  is  a  stability  theorem  for  a  certain  class  of  dif¬ 
ference  algorithms  designed  to  give  approximate  solutions  of  a  model  inverse  scattering 
problem  in  one  dimension.  This  stability  result  guarantees  the  convergence  of  the 
approximate  solutions  to  the  exact  solution  of  the  problem  as  the  grid  of  the  difference 
scheme  is  refined.  We  present  the  results  of  numerical  experiments  based  on  one  of 
these  schemes,  in  which  second-order  convergence  is  observed.  Furthermore  the  cost 
(that  is,  the  dependence  on  N  of  the  number  of  arithmetic  operations  required  to 
compute  the  solution  at  N  grid  points)  of  the  algorithms  discussed  below  is  essen¬ 
tially  optimal. 

The  algorithms  of  this  paper  are  difference  approximations  to  a  certain  hyperbolic 
boundary  value  problem.  In  a  previous  paper  |1),  we  showed  that  this  hyperbolic 
boundary  value  problem,  the  solution  of  which  leads  to  a  solution  of  the  above-mentioned 
inverse  scattering  problem,  is  equivalent  to  a  certain  Volterra  integral  equation,  and 
the  latter  was  solved  constructively,  with  estimates.  Roughly  the  same  plan  is  followed 
in  this  paper.  An  approximate  version  of  the  Volterra  equation  is  first  solved,  with 
estimates.  This  discrete  Volterra  equation  is  then  shown  to  be  almost  equivalent  to  a 
certain  difference  approximation  to  the  hyperbolic  boundary  value  problem.  The  relation 
is  close  enough  that  stability  estimates  for  a  class  of  difference  approximations 
follow. 

It  is  interesting  that  the  usual  approaches  to  proving  stability  of  difference 
schemes  for  time-dependent  problems  (see  for  instance  (2))  do  not  seem  to  apply  to  the 
problem  considered  here. 

Sponsored  by  the  United  States  Army  under  Contract  No.  DAAG2'9-'K-C-0024.  This  material 
is  based  upon  work  supported  by  the  National  Science  Foundation  under  Grant  No. 
MCS78-09525 . 


$2.  Statement  and  Discussion  of  Results 

We  fust  discuss  the  continuum  problem  and  its  solution,  with  various  estimates, 
as  presented  in  11).  The  continuum  results  provide  the  major  tools  for  the  approximate 
solution,  as  well  as  the  form  of  the  estimates  which  must  be  reflected  in  the 
stability  theorem  for  the  various  difference  schemes.  We  then  present  some  notational 
conventions  (concerning  difference  schemes  and  uniform  estimates)  to  which  we  will 
adhere  throughout.  Finally,  we  state  the  difference  schemes  which  form  the  main 
subject  of  the  paper,  and  outline  their  stability  properties. 

The  inverse  scattering  problem  of  this  paper  is:  given  a  real-valued  function 
F  i  (0,  2T J  »,  find  a  real-valued  V  s  (0,  T)  v  H  and  H  c  1R  so  that 
F ( t )  *  U(0,  t) ,  0  <  t  2T,  for  the  solution  U  of  the  initial  boundary-value 
problem 

/  a2  a2  \ 

— -j - X  *  V (x) /  U(x,  t)  -  0  (2.1a) 

\3t  3x  / 

(lx  +  HU)  (0<  "  0  t  >  0  (2.1b) 

U(x,  0)  -  C$(x),  —  (X,  0)  S  0  (2.1c) 

In  the  region  {(x,  t):  0<x<I,  0<^t<  2T).  We  shall  refer  to  this  problem  as 

the  inverse  problem. 

This  inverse  problem  can  be  considered  a  very  simple  instance  of  the  inverse 
scattering  problem  for  a  mechanical  continuum  supporting  small-amplitude  wave  propaga¬ 
tion.  Roughly  speaking,  you  are  required  to  construct  a  vibrating  one-dimensional 
continuum  having  equations  of  motion  of  the  form  (2.1a),  boundary  condition  of  the 
form  (2.1b),  and  having  a  prescribed  response  (back-scattered  wave  at  x  ■  0) 

F(t),  t  »  0,  to  an  impulsive  ("broad  band")  incident  disturbance  (initial  conditions 
(2.1c)),  for  a  prescribed  duration  2T  >  0.  The  geometry  of  such  a  hyperbolic  mixed 
problem  shows  clearly  that  the  boundary  value  F(t)  -  11(0,  t)  for  0  <  t  2T 
depends  on  the  coefficient  V(x)  only  for  0  x  T. 
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Thi»  problem  is  typical  tn  various  re*i>ectx  of  a  large  mimN-i  of  time-dependent 
aiuI  time- independent  inverse  scattering  problems. 

Various  other  problems  which  may  be  treated  by  the  techniques  descr  ibe.t  below  ai> 
discussed  briefly  in  the  Conclusion  (section  tt) .  furthermore,  at  least  tn  the  limit 
r  *  "•  this  problem  has  a  long  and  illustrious  history,  bsinq  a  time-dependent 
version  of  the  inverse  spectral  problem  for  Sturm-Liouville  operators  solved  by 
I.  M.  Gel’fand  and  H.  M.  Levitan  in  their  seminal  paper  11).  The  connection  t<etween 
the  half-axis  version  of  the  inverse  problem  for  (la,b,c)  and  the  Gel ' fand-Levitan 
inverse  spectral  problem  is  explained  in  11),  section  3, 

It  was  shown  in  11)  that  solution  of  the  inverse  problem  is  equivalent  to  solu¬ 
tion  of  a  hyi>erbolic  boundary  value  problem  1(2.2)  below)  .  The  Russian  mathamat ioian 
Chudov  seems  to  have  been  the  first  to  suggest  the  possible  use  of  this  boundary 
value  problem  as  a  means  for  solvinq  inverse  problems  (see  (2),  final  section).  We 
shall  therefore  call  this  problem  the  Chudov  problem.  As  we  shall  explain  elsewhere, 
the  approach  to  inverse  problems  through  the  Chudov  problem  is  closely  related  to  the 
recant  work  of  Delft  and  Trubowit*  15)  on  the  1-dimeneional  Sohrbdlnger  inverse 
scattering  problem,  and  to  the  wor*  of  Hochatadt,  Hald,  and  Levitan  (|t>|,  17),  )H))  on 
inverse  Sturm-Liouvil le  problems. 

The  Chudov  problem  is  as  follows.  Let  C^,  -  {  (x,  t)  i  0  <  x  <_  min(t,  2T  -  t)). 

find  W  i  C_  ■»  K  such  that 
T 

/  \ 

(  — 3  -  — T  ♦  V(x)  I  N(x,  t)  -  0  (2. 2a) 

Vat  a« 
w(o,  t)  -  r(t) 

1 5^  ♦  HW )  (0,  t)  -  0  2T  >  t  »  0  (3.2b) 

h  -  r(0) 

Vlx)  -  -  2  -  W (X,  x)  0  <  x  <  T  (2.2c) 
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Here ,  as  In  the  Inverts  problem,  the  datum  of  the  problem  1»  the  function 
r  i  10,  2T  1  -*  *.  Note  that  12.2)  ts  nonlinear,  by  virtue  of  the  coupling  of  solution 
ami  coefficient  in  (2.2c). 

In  our  previous  sank  (II,  estimates  are  g i ve n  for  various  Sobolev  mums  of  N, 

V  m  terms  of  similar  norms  of  F,  and  a  number  i  '  0  with  the  property  that,  for 
any  ♦  •  L*  (0,  T), 

T  3  T  T  .  T 

/  |#|*  ♦  /  ds  /  dt  ♦  (*)  »(t)  <  J  (F  (a  ♦  t)  ♦  r  1 1  •  -  1 1 » >  >  .  /  >l!  (3.1) 

0  0  0  0 

The  results  of  |1)  ate  phrased  in  terms  of  the  spaces  N* ' ^  of  functions  with 

n  absolutely  integrable  derivatives.  Another  somewhat  unusual  Sobolev  Space, 

denoted  by  tt2*.l 

in  111,  intervenes  as  an  auxiliary  device.  This  space  ia,  roughly 
speaking,  the  subst'ace  of  w** ’  i  11* )  in  which  each  function  has  a  well-defined 
restriction  to  each  line  in  K^,  lying  in  tfB|  *  of  functions  on  the  line,  and  the 
restrictions  vary  continuously  with  the  choice  of  line.  Precise  definitions  ate  found 
in  11),  section  2. 

The  following  result  is  proved  in  sect  tons  5,  *  of  ( 1 ) t 

2  1 

Theorem  1,  Problem  (2a,b,c)  has  a  solution  N  in  the  space  It'  ’  (C_)  if  and  only 

if  F  satisfies 

(i)  F  <  NJ*1  (10,  2X1) 

(ii)  there  exists  •  >  0  so  that  (2.))  is  satisfied  for  all  ♦  <  l.”  (10,  T )  > . 

If  these  are  satisfied,  then  V  «  W*’1  1(0,  Tl)  and  ia  given  by  (2.2c). 

It  is  convenient  to  paraphi aae  condition  (ii)  of  the  theorem  as  follows.  Define 
the  symmetric  kernel 

Gte,  t)  -  '  (F<a  ♦  t)  ♦  F 1 1 s  -  t|)> 

If  F  satisfies  (i) ,  the  kernel  0  ia  continuous,  hence  defines  a  bounded  self- 
adjoint  Nilbert-Schntidt  operator  G  on  L*  (10,  Tl).  Denote  by  1  the  identity 
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i-nliWlKV  TiiliWftW 


operator  on  L*  (  |0,  T|).  Then  (ii)  can  be  stated: 

I  ♦  G  >  t  >  0 


(ii*) 


The  method  of  solution  of  the  Chudov  problem  (2.2)  is  based  on  the  following 
result : 

Theorem  2.  A  function  W  e  li' ' 1  (C  )  solves  the  Chudov  problem  (2)  if  and  only  if 
-  T 


Gls ,  t)  *  W(s,  t)  ♦  /  dy  W(y,  s)  w(y,  t)  for  (s,  t)  t  C  . 


(2.4) 


Remark  1)  We  shall  refer  to  (2.4)  as  the  G-L  equation,  or  (GL)  .  It  appears  first 
in  the  paper  of  Gel'fand  and  Levitan  13),  and  also  expresses  the  group  law  of 
propagation  of  initial  values  for  the  mixed  problem  (la,b). 

2)  for  (s,  t)  t  C  ,  we  have  0  <  s  <  t.  If  we  associate  to  the  kernel  W 

T  —  — 

(which  is  at  least  continuous)  the  Volterra  operator 


W  ♦(y)  -  /  dt  W (y ,  t)  ♦  (t) 

y 

then  you  can  write  the  GL  equation  in  the  form 

I  ♦  G  -  (I  4  W)  +  (I  ♦  W) 

which  clearly  illustrates  the  necessity  of  hypothesis  (ii*)  of  Theorem  1. 

The  GL  equation  allows  the  derivation  of  a  number  of  estimates  for  V  in  terms 
of  F  in  various  norms  other  than  the  Sobolev  (m,  1)  -norms  of  the  (sharp)  stability 
statement.  For  our  purposes  the  following  C^  -estimates  will  be  sufficient. 


Ilwll  <  ft  Fll  (1  +  t’1  T  IlFll  ) 
®  —  ®  » 


II D  Wll  <  IlDFll  exp  (T  II  Wll  ) 

q  »  —  * 

IVll  <2  IlDFll  +  IIWlI  (Ilwll  +  2T  II  D„  Wll  ) 

OB  —  tv  (V  v»  J  vV 


(2.5a) 

(2.5b) 

(2.5c) 


At  this  point,  it  should  be  mentioned  that  the  estimate  (2.5)  is  not  really  a 
stability  result,  since  the  problem  (2.2)  is  nonlinear.  However  it  is  very  easy  to 
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prow  genuine  stability  results  (that  is,  local  Lipschitz  estimates  for  the  map 
F  **  V,  H)  on  the  basis  of  boundedness  results  like  (2.5).  For  both  the  continuum  and 
approximate  problems,  therefore,  vie  will  refrain  from  stability  statements  (which  are 
necessary,  for  instance,  to  derive  explicit  error  bounds)  and  present  only  boundedness 
results  like  (2.5). 

This  concludes  our  discussion  of  the  (continuous)  inverse  problem.  We  now  turn 
to  approximate  methods  of  solution,  and  begin  by  establishing  notation  and  terminology’ 
for  the  difference  schemes  we  shall  use. 

The  approximate  algorithms  of  this  paper  are  difference  schemes  for  computing 
certain  grid  functions.  The  grids  are  uniform  and  rectangular,  and  the  granularity 
will  be  denoted  by  A.  We  will  use  the  same  letters  for  the  discrete  approximants  as 
for  their  continuous  counterparts;  thus,  F(n)  is  meant  to  approximate  Find), 

W(n,  m)  corresponds  to  W(nA,  mA),  etc.  The  number  of  (linear)  gridpoints  will  be  N. 


We  will  have  need  of  the  sup  norms 


II  Ftl  «  sup  |F(n)| 
l<n<N 

II  Ml  «  sup  |w(n,  m)  | 
l<n,m<N 


and  occasionally  the  norms 


n  1/p 

DFtl  ^  -  (  l  A  |F(u)  |p) 
p  n-1 


1  <  p  < 


The  basic  difference  operators  defined  on  grid  functions  are  given  by 
D+  F(n)  -  A_1(F(n  *  1)  -  F (n) ) 

D*  F(n)  -  A-1(F(n)  -  F (n  -  1) ) 

The  partial  difference  operators  on  2-dimensional  grid  functions  will  be  denoted  by 
subscripts,  for  instance: 


Dj  W(n,  m)  -  A  Swin  +  1,  m)  -  W(n,  m)) 
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and  so  on.  Wo  will  also  need  the  diagonal  derivatives 


cXir  aim  la  to  produce  estimates  for  various  difference  schemes,  based  on  such  grids 


which  are  uniform  in  A  as  A  »  0,  that  is,  with  T  •*  NA  fixed 


The  first  set  of  estimates  i>ertain  to  solutions  of  an  approximate  GL  equation 


!T)  -»  K  is  a  Lipschit*  function,  and  set 


We  sup<posv  that  F 


If  we  approximate  the  integral  in  (GL)  by  right-endpoint  Riemann  sums,  we  obtain  the 


discrete  Volterra  equation 


for  1  <  n  <  m.  We  can  require  (harmlessly) 


In  section  4,  we  obtain  the  following  results: 


1)  The  system  (2.6  -  2.7)  can  be  written 


where  I  is  the  N  *  N  identity  matrix,  and  thus  represents  the  Cholesky  decomposition 


I  >  G.  A  solution  therefore  exists  if  and  only  if  the  L.H.S.  is 


of  the  matrix 


positive-def inite,  with  lower  bound  t  >  0 


2)  Under  suitable  conditions,  which  we  shall  not  discuss  here,  the  solution  w  of 


(2.6)  -  (2.7)  converges  in  sup  norm  to  the  solution  of  (2.4) 


3)  The  solution  W  of  (2.6)  -  (2.7)  is  estimated  by  (Proposition  4.1) 

UUI  <  HF«  ♦  t'1  ■  Ffl  *  (2.8) 

Further,  provided  4  is  small  enough  in  relation  to  II  Fll  .  one  can  also  bound  the 
partial  differences  of  W  by  entire  functions  of  IIWII  and  norms  of  differences  of 
F,  which  are  linear  in  F  near  F  •  0  (Propositions  4.2,  4.3). 

It  is  evident  from  (2.8)  that  uniform  estimates  as  the  grid  is  refined  (4  ♦  0) 
will  only  be  achieved  if  the  lower  bound  t  •  t (4)  eventually  stabilises.  Our 
results  therefore  apply  when  the  discrete  approximations  to  the  backscattered  wave 
have  been  chosen  so  as  to  guarantee  such  a  uniform  lower  bound.  We  will  not  discuss 
methods  for  extracting  such  discrete  grid  functions  from  experimental  data  in  this 
paper,  although  this  very  interesting  matter  should  certainly  be  considered  further. 

Various  elementary  estimates  and  identities,  used  to  prove  the  above-mentioned 
estimates  and  for  various  other  purposes,  are  collected  in  section  3  for  easy 
reference. 

The  identities  and  estimates  of  section  3  are  used  in  section  5  to  compute  a 
difference  scheme  satisfied  by  the  solution  W  of  (2.6)  -  (2.7).  The  result  is 
(equation  (5.14)): 


(D*  D”  -  D*  D~)W  ♦  VW  -  4R 


(2.9a) 


V  (n) 


(o'  +  D*  ♦  D~)W(n,  n) 


(2.9b) 


D*  W(0,  m)  ♦  HW(0,  m)  »  4B(m,  A) 

W(0,  m)  ■  F(m) 

Here  R  (2-dimensional)  and  B  (1-dimensional)  are  grid  functions  whose  sup  norms 
are  bounded  in  terms  of  II  Fll  ,  II D-  Fit  ,  IlD*  D  Fll  ,  T,  4,  and  t. 

This  result  is  useless  for  computation  purposes,  since  the  R.H.S.  of  (2.9a), 
(2.9c)  depend  explicitly  on  W  in  a  conflicated  way.  On  the  other  hand,  (2.9)  is 
clearly  related  to  the  following  difference  approximation  to  (2.2): 


} 


(2.9c) 
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(D2  °2  '  °1  °1)W0  *  V0  W0  “  ° 


(2.10a) 


V  (n)  -  -  (D  >  D  +  D,  )W  (n,  n) 
0  n  2  1  0 


(2.10b) 


W^(0,  m)  «  F(m) 

D+  WQ<0,  m)  ♦  H  W^IO,  m)  «  0 


(2.10c) 


The  mam  result  of  this  paper  is  the  following  estimate,  proved  in  section  7; 
relating  the  solutions  of  the  systems  (2.9)  and  (2.10): 

II W  -  W0II  <  A  C 

(2.11) 

II  v  -  vj  <  a  c„ 

o  -  2 

Here  C^,  are  entire  functions  of  II  Fll  ,  II D—  Fll  ,  II D+  D  Fll  ,  T,  «  and  A.  (In 
fact,  C^,  are  exponential  polynomials;  and  can  be  written  out  explicitly, 

although  ve  do  not  do  so  here). 

The  estimates  (2.8)  and  (2.11)  combine  to  yield  boundedness  statements  like 
(2.5)  for  the  solution  of  (2.10).  An  interesting  difference  is  that  the  resulting 
estimates  on  II VHI  involve,  in  the  limit  A  -*  0,  the  modulus  of  Lipschitz  continuity 
of  the  derivative  DF,  whereas  the  sup  norm  estimate  for  V  in  the  continuum  problem 
(2.5c)  only  requires  that  DF  be  continuous. 

The  system  (2.10)  is  still  inefficient  for  computation;  as  it  has  local  trunca¬ 
tion  error  (on  the  boundaries)  of  first  order.  A  second-order-consistent  approxima¬ 
tion  to  (2)  is,  for  instance. 


(D*  D~  -  D*  D~)W  +  VW  -  0 


(2.12a) 


(0,  m)  *»  F  (m) 

W  (1,  m)  =  (jV(0)A2  -  HA) W (0 ,  m)  +  i  (W(0,  m  +  1)  +  W(0,  m  -  1) ) 


(2.12b) 


V(n)  »  -  (D  +  D  )W(n,  n) 
n  n 


(2.12c) 
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A  bonus  of  our  method  is  that  the  stability  of  the  lst-order-accurate  system  (2.10) 
implies  the  stability  of  the  2nd-order-accurate  system  (2.12).  Indeed,  the  two  differ 
by  terms  which  are  0(A),  which  may  be  shifted  to  the  R.H.S.  Our  arguments  are 
exactly  designed  to  extend  the  estimates  for  (2.10)  to  similar  estimates  for  (2.12). 

Well-known  arguments  then  guarantee  that  the  discrete  coefficient  V,  constructed 

by  solving  (2.12)  on  a  machine,  differs  from  the  exact  solution  of  (2.2),  evaluated 

2 

at  the  corresponding  grid  points,  by  an  error  proportional  to  A  .  The  constant  of 
proportionality  can  even  be  estimated  in  terms  of  sup  norms  of  differences  of  F,  T, 
e  ^ ,  and  the  round-off  characteristics  of  particular  machines,  though  we  shall  not 
do  this  here. 

Instead,  we  give  the  results  of  some  numerical  experiments  based  on  the  second- 
order  scheme  (2.12).  These  are  displayed  in  section  7.  We  end  with  a  discussion,  in 
section  8,  of  related  problems  which  may  be  solved  by  the  same  methods,  as  well  as  the 
relation  of  our  results  to  previous  work  on  inverse  scattering  problems. 
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1'\rcr  t  "..'V*  %  .  . 


i).  Elementary  Identities  and  Estimate* 

Wo  begin  with  soino  estimates  for  discrete  Vo  L  ten  a  pqvwt  Uhib  (obtained  ftom 
Vol terra  aquat  ions  by  raplacing  integrals  with  right  -eint*point  Kiomami  sums)  .  These 
aro  entirely  elementary,  but  wo  know  of  no  reference  which  ntaton  exactly  the  esti¬ 
mates  wo  need,  so  wo  give  complete  proofs.  Wo  end  with  a  number  of  useful  identities 
from  the  calculus  of  finite  differences, 
bomma  (J.l)  Suppose  that,  for  1  <  n  «.  N 


gin)  -  ♦  (n)  *  J  i  W(n,  k)*(k) 
k-1 


Mil  <  Hull  ext'  <NA  II  Mil  ) 


Proof.  Define  TV (n)  *  g(n)  -  £  A  W(n,  k)\Mn)  for  any  grid  function  iiMn), 

k-1 

1  f.  n  i.  N)*  Then  wo  seek  a  Tixed  point  of  T.  The  point  is.  of  course,  that  the 
fixed  i<oint  exists  and  is  globally  asymptotically  stable.  For  instance,  sot 

♦o  *  ° 

Vi " T  V 

and  u  •  $  *  ♦ 

p  rp  Tp-i 

Tht»n  u  .  •  T.  u  ,  whf>re 

p>  1  Op 


T0  g. in)  -  -  l  A  win,  kX'(k) 

U*1 


(nA  llWtl  J1 

lup(n)  I  - - pT*~~  llg*w  p  •  l*  2>  *  *  * 

The  claim  certainly  holds  for  p  -  1.  Suppose  that  it  holds  up  t o  p  -  1.  Then 


n-r 

-  T  A  W(n,  k)u  ,  Ik)  | 

k-1 


11 


n-1 

g(n)  -  ^(n)  ♦  £  A  Win,  K ) <|> ( k ) 

k-1 

hence  according  to  the  previous  lemma, 

II  i|HI  ^  <_  II  gll  exp  (NA  II  Wll  m ) 

which  immediately  implies  the  asserted  estimate. 

q.e.d. 

Lemma  Suppose  |U(n,  m)  :  l<^n<^m<^N}  satisfies 

n-1  n-1 

G(n,  m)  »  U(n,  m)  +  £  A  W  (k,  n)l)(k,  m)  +  J  A  U(k,  n)W  (k,  m) 

n-1  1  k-1  1 

Then 

Hull  II  dl  exp  |NA  max  (llw^ll,  II  W  Jl  )  ] 

Proof.  Exactly  parallel  to  the  proof  of  Lemma  1. 

q.e.d. 

We  conclude  with  a  number  of  "summation  by  parts"  formulas  which  will  be  used  in 
section  5.  The  proofs  are  all  trivial,  so  we  omit  them. 

1.  D+  d"  -  D_  D+ 

n-1  n-1 

2.  I  A (D  f(k))g(k)  -  -  £  A  f (k)D~  g(k)  +  f(n)g(n  -  1)  -  f(l)g(0) 

k-1  k-1 

n-1  n-1 

3.  I  aid'  f(n))g(k)  -  -  l  A  f(k)  (D  g(k>)  +  f(n  -  l)g(n)  -  f(0)g(l) 

k-1  k-1 

4.  D"  f  (k)  -  Df  f  (k  -  1)  ,  D_  f  (k  +  1)  -  D+  f(k> 

n-1  n-1 

5.  I  A(D  D  f  (k) )  g  (k)  -  £  A  f(k)  ID  D  g(k)) 

k-1  k-1 

+  D_  f(n)gln)  -  f(n)D_  gin)  -  D_  f(l)g(l)  +  f |1)d"  g(l) 
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S4.  Estimates  for  Cholesky  Decompositions 

bet  lF(n).  n  »  0,  •  *  •  ,  N)  Iw  an  (N  ♦  1)  -  vector,  and  define  the  symmetric 
array  G  by 

G(n,  in)  •  ^  (F(n  +  m)  4  F(|n  -  m|)J 
1  -»N 

Suppose  that  the  matrix  (  -r  I  +  G|  ,  is  positive-definite,  where  I  is  tlie  N  »  N 

t  A  Jn.iiel 

identity  matrix.  Then  ^  X  +  G  admits  a  Cholesky  decomposition 

(}  i  +  c)  •  a  4  i  ♦  wS  (|  i  ♦  w)  (4.x) 

A  A  A 

where  W  is  a  triangular  array,  W(n,  m)  •  0  if  n  >  m.  Xn  this  section  we  obtain 

estimates  for  w  and  its  partial  differences  in  terms  of  F  and  its  differences,  and 

the  lower  bound  for  7  I  i  C, 

A 

The  formula  (4.1)  can  1-e  rewritten 

n 

G(n,  m)  «  W(n,  m)  4  £  A  W(k,  n)W(k,  m)  (4.2) 

k-1 


for  m  >  n.  This  suggests  the  introduction  of  the  Hilbert  subspaces  H  «  l(  • 

m 

t‘ll,  •  •  •  ,  N)  defined  by 
A 

g>  c  H  •  (Mn)  “0  ,  n  >  m 

m 


bet  n  i  H  ■+  H  t'e  the  orthogonal  prelections.  Then  (4.2)  can  l->e  written 
m  m 

n  g  -  n  a  (t  +  w+ )  n  w 

m  n  m  A  mm 

(here  G  (n)  »  G(n,  m) ,  W  (n)  «  W(n,  m) ) 
m  m 

From  (4.1)  you  obtain 


11  n  c  11  -  <w  ,  n  a  (i  +  win  (i  +  wSn  w  > 

mm.  m  m  A  m  A  mm 


(4.1) 


(4.4) 


where  <  ,  >  is  the  inner  product  associated  with  the  norm  II  II  . 

We  claim  that  the  operator  in  the  R.H.S.  of  (4.4)  is  Invertible.  Xn  fact,  it  is 
of  the  form  K  K  ,  where 
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k  -  n  (7  +  win 

m  A  m 

is  an  operator  on  H  .  Notice  that  the  triangularity  of  W  implies  that 
n\ 

n  w  n  -  w  n 

mm  m 


k+  k  -  a  n  (7  +  w'  in  n  <7  +  w>n 

tn  A  m  tn  A  m 

-an  (7  +  n+  i  (7  +  w)nm 
m  A  A  m 

-  n  (7  +  Gin 
m  A  m 

Since  —  +  G  is  bounded  below  on  H  by  e  >  0,  it  is  a  fortiori  bounded  below  by 

t  when  sandwiched  between  projections  on  H  .  It  follows  that 

m 

II  K  *11  2  ^  t  II  *11  2 

for  all  *  t  H  .  Then  K  must  also  be  invertible,  hence  K  K  is  invertible, 
in 

Since  both  K  K  and  K  K+  are  invertible  on  H  ,  they  have  the  same  spectrum, 

m 


which  necessarily  lies  above  e .  Thus 


(w  ,  n  a  (7  +  win  (7  +  w+m  w  ) 

mm  A  m  A  m  in 

<n  w  ,  tn  a  (7  +  win  (7  +  w+nn  w  1 

mmmA  mA  mm 


>  t  •  lln  W  H  ‘ 
—  mm. 


If  n  <  m,  then  w  -  n  W  and 
n  m  n 


(W  ,  W  )  -  l  A  W(k,  n)W(k,  ml 


<n  w.n  w  > «  <n  w,n,wi 

tn  n  m  tn  n  n  tntn 


On  the  other  hand 


n  n  c  ir’  <  ii  m  * 

mm,—  2 


(4.7) 


By  conbining  (4.4),  (4.5),  I4.t>),  and  (4.7)  and  using  the  Cauchy-Schwartz  inequality 
one  obtains  finally  that  for  n  <  m 


(W  .  W  >  <  f'1  II  Fll  ^ 
n  m  —  2 


(4.8) 


and  so  from  (4.2)  that  for  n  <  m 


|w(n,  m)  |  |G(n,  m)  |  +  t  1  II  Fll  * 


(4.8) 


To  estimate  the  diagonal  elements  W(n,  n) ,  write  from  (4.1) 

n-1 


G(n,  n)  ■  2W(n,  n)  +  A  W(n,  n)*  +  \  A  W(k,  n)' 

k-0 


Thus 


1  /  "'1  2  1 
w(n,  n)  »  +  —  »  1  +  lG(n,  n)  -  7  A  W(k,  n)  )A  -  — 

—  A  **  A 


k«0 


In  order  to  maintain  the  positivity  of  the  diagonal  elements  of  —  1  +  W,  as  is 

A 


required  of  the  Cholesky  decomposition,  we  must  choose  the  +  sign  in  front  of  the 
radical.  Since 

n-1 


|/l  +  (G(n, 


n-1  ,  |  |  n-i  | 

n)  -  l  A  W(k,  n)  )A  -  l|  <  |G(n,  n)  -  J  A  W(k,  n)  |  A 

k-*  "  k-t 


Another  application  of  (4.8)  shows  that  (4.9)  is  valid  also  for  n  ■  m.  So  one  has 
proved 

Proposition  4.1.  The  cholesky  factor  W  satisfies 


IIWII  <  It  Fll  ♦  t'1  II  Fit  \  <  « Fll  (1  +  t*1  N  A  IlFll  ) 
®  —  00  2  «  ® 


The  next  step  is  to  estimate  the  partial  differences  of  W.  We  start  with 


W(n,  m)  «  i  lW(n,  m  +  1)  -  W(n,  m)  I  ,  n  <_  m 
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Now 


which,  together  with  Propositions  4.1  and  4.2,  gives  (a)  +  .  Similar  computations 
give  the  other  estimates. 
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15.  Partial  Difference  Equation  for  Cholesky  Decomposition 


W  will  denote  the  lower  triangular  solution  to 


As  in  the  previous  section 


the  first  summand  on  the  R.H.S.  of  (5.1)  can  be  dropped,  thus 


For  the  remainder  of  the  paper,  we  use  the  symbol  0  to  abbreviate  the  usual 


five-point  approximation  to  the  wave  operator 


An  easy  calculation  shows  that 


Therefore,  applying  0  to  both  sides  of  (5.3)  and  using  formula  6  (from  section  3) 


gives 


n- 1  ,  *»  *  + 

0  M(n.  m)  ♦  l  A  W(k,  n)  (D*  D~  W(k.  m) )  -  V  A(D  D.^  W(k,  n))W(k,  n\) 


ID*  W(n,  n)  ♦  D*  Win,  n))W(n,  m)  »  -A  OlWln,  n)W(n,  m) )  ♦  Win  -  1,  n  -  1)D1  Win,  m) 


Also  formula  5  of  section  3  implies 


n-1  .  n-i  _  + 

y  AID*  D,  W(k,  n))W(k,  m)  -  5  A  W(k,  n)  (D,  D,  W(k,  m) )  +  ID  W(n,  n))Wln,  m) 
,  ,  1  1  .  ,  11  1 


-  Win,  n)D~  W(n,  m)  -  (D ^  Wll.  n))W(l,  m)  +  Wll,  n)Dj  Wll.  m) 


We  shall  examine  the  last  two  summands  in  the  above  more  closely.  For  n  *  1  <  m. 


(5.1)  reads 


G(1 ,  m)  -  ~  (F (m  +  1)  +  Flm  -  1))  »  (1  ♦  A  W(l,  1))W(1,  m) 


and  for  n  -  m  «  1 


i  IF  (2)  +  F  (0)  )  -  W(l,  1)  +  A  Wll,  1) 


It  follows  that 


W(l,  1)  -  F (0)  +  0^  l A )  -  H  +  (A) 


where  the  bound  x  in  the  0.  statement  can  be  estimated  in  terms  of  I!  Fll  ,  IlDFtl  . 

1 


On  the  other  hand. 


j  IF (m  +  1)  +  F(m  -  D)  -  F(m>  +  ~  D+  d'  Flm) 


Combine  these  facts  to  obtain 


Wll,  m)  -  (1  -  HA)F(m)  +  0.,(A) 


where  the  bound  in  0^  can  be  estimated  in  terms  of  IlFll^,  IlDFll^,  II D  D  Fllw.  It 


follows  (recalling  that  W(0,  m)  *  F(m))  that 


D  Wjll,  m)  -  -H  Flm)  +  C^IA) 


and  finally  that 


1 


(o'  Nil,  n)  )W(1,  m)  -  W(l,  n)Dj  W(l,  m)  -  -H  F(n)  (1  -  HA) F (in) 

♦  H  F(m)(l  -  HA)F(n)  +  Ox  (A)  5  -A  Rj(n,  m.  A) 

where,  for  A  sufficiently  small,  llR.II  can  be  estimated  in  terms  of  II  Ftl  , 

1  •» 

IlDFll  ,  and  II D+  D~  Ftl  . 

•  * 

Now  rewrite  (5.6)  as 

n- 1  n-1 

l  A(D'  D*  W (k ,  n)  )W(k,  m)  -  l  A  W(k,  n)  <D~  D*  W(k,  m) ) 
k-1  11  k-1 

(5.9) 

-  (D^  W(n,  n))W(n,  m)  -  -W(n,  n)D^  W(n,  m)  +  A  R^(n,  m,  A) 

Add  this  identity  to  (5.5)  to  obtain 

n-1  n-1 

0  W(n,  m)  +  £  A  W(k,  n)0  W(k,  m)  -  £  A  (0  W(k,  n))W(k,  m) 

k-1  k-1 

(5.10) 

♦  V(n)W(n,  m)  =  A  R(n,  m,  A) 

V(n)  -  -(D*  +  D*  +  D")W(n,  n) 
nil 

where 

R <n »  *  -0(W(n,  n)W(n,  m)J  -  Dn  W(n,  n^  W(n,  m)  +  R^n,  m,  A)  (5.11) 

According  to  the  results  of  section  4,  II  Rll  may  be  estimated,  for  small  A,  in 

terms  of  II  Fll  ,  IlDFll  ,  II D+  D~  Ftl  . 

C®  CP  OD 

Notice  that  we  can  replace  the  factor  0  W(k,  m)  in  the  first  sum  on  the  L.H.S. 
of  (5.10)  by  0  W(k,  m)  +  V(k)W(k,  m) ,  provided  we  also  replace  0  w(k,  n)  in  the 
second  sum  by  0  w(k,  n)  +  V(k)W(k,  n) . 

Define 

U(n,  m)  -  0  W(n,  m)  +  V(n)W(n,  m)  ,  1  n  <  m 

Then  U  obeys  the  discrete  Volterra  equation 
n-1  n-1 

U(n,  m)  +  £  A  W(k,  n)U(k,  m)  -  £  A  U(k,  n)W(k,  m)  =  A  R(n,  m.  A)  (5.12) 

k-1  k-1 

According  to  Lemma  3.3,  the  solution  enjoys  the  same  sort  of  bounds  as  the 
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inhomogeneous  term.  In  particular,  we  get  1)  *  AP,  so  that 


C1  W(n,  m)  +  V(n)W(n,  m)  =  A  P(n,  m.  A)  (5.13) 

where  II  Pll  m  is  bounded,  uniformly  for  small  A,  in  terms  of  IlFll^,  II  DEll  ^ ,  II D+  D~  Fll  , 
and  l|W||  ,  hence  (according  to  the  results  of  section  4)  in  terms  of  II  Fll  ,  II  DFll  , 

00  OB 

II D  D  Ftt  and  e 

For  future  reference,  we  collect  (5.13),  (5.11),  and  (5.8)  in  the  form  of  a 
boundary- value  problem  for  W: 

0  W  +  VW  =  A  P  (5.14a) 


V(n)  =  -(D  +  D  +  D  ) W (n ,  n) 

n  2  1 


W(0,  m)  =  F (m) 


Dx  W(0,  m)  +  H  W(0,  m)  =  A  B(m,  A) 


(5.14b) 


(5.14c) 


where  in  the  last  line  of  (5.14c),  the  vector  B(m,  A),  M  >_  0,  defined  implicitly  in 
the  foregoing,  is  bounded  uniformly  for  small  A  in  terms  of  II  Fll  ,  II  DFll  ,  and 


§6.  Estimates  for  Chudov  Schemes 

The  purpose  of  the  present  section  is  to  compare  the  solution  W  of  the  Cholesky 
equation  (5.1),  (alternatively,  of  the  boundary-value  problem  (5.14)),  to  the  solution 
Wq  of  the  first  order  Chudov  Scheme  (equations  (2.10)  of  section  2),  which  we  display 


here  for  convenience: 


0  W  +  V„  W  =  0 
0  0  0 


VQ(n)  =  -<Dn  +  D2  +  Dx)Wo(n,  n) 


WQ(0,  m)  =  F(m) 


D  WQ(0,  m)  +  H  WQ(0,  m)  =  0 
H  =  F  (0) 

The  following  notational  convention  is  convenient:  for  any  gridfunction 

A(n,  m)  defined  for  n  <_  m,  set 

D  A (n,  n)  =  -  (D~  +  D*  +  D')A(n,  n) 
n  2  1 

Then  (6.1a),  (6.1b)  together  become 

0  WQ(n,  m)  +  WQ(n,  m)D  WQ(n,  n)  «  0 

Similarly  (5.14a),  (5.14b)  become 

0  W(n,  m)  +  W(n,  m)D  W(n,  n)  =  A  P(n,  m.  A) 

Subtracting  (6.2)  from  (6.3)  we  obtain  for  the  difference  W  =  W  -  W 


0  w(n,  m)  +  V(n)W(n,  n)  +  w(n,  m)D  W(n,  n)  -  W(n,  m)5  W(n,  n)  =  A  P(n,  m,  A)  (6.4) 
W  also  satisfies  the  boundary  equations 

W(0,  m)  =  0 


Dx  W(0,  m)  +  H  W(0,  m)  =  A  B(m,  A) 

We  shall  show  that  both  W  and  its  first  differences  are  0(A)  as  A  -*•  0,  in  the 


where  N  =  T/A. 


CT  =  {  (n,  m)  :  0  £  n  <_  min  (m,  2N  -  m)  } 


,  •  .... _ 


_ _ _ _  _ , _ . 


Til*  main  toot  Is  the  next  lemma,  which  show*  that  hounds  of  t  h<>  type  wv  want 
lx>  piopauat  ed  to  finite  "depth"  utA  const . )  ,  with  controlled  growth  tn  the 
constants  of  proport ional ity. 

l.emma  (h.l)  suppose  that  W  satisfies  the  difference  equation 

0  W(n.  m)  ♦  V(ttlW(n,  ml  twin,  mil'  Win.  n)  -  Wtn,  mlf1  Witt,  n\  -  A  Pin,  m,  A\  , 

m  the  half-lattice  (lit,  m)  t  n  '  0}.  Suppose  tn  addition  that  the  fol  lowing 

estimates  hold: 

|W(0,  m) |  ,  |W(1,  m) |  '  A  C 

| W(  1 ,  m)  -  W(0,  m  -  1 )  |  v  A^  t'( 

for  m  «  X. ,  aihI 

|w<n,  ml  |  <  Cj  ,  | V in)  |  • 

| Pin,  m)  |  ^  C 

for  n  ^  0,  m  «  *.  Then  for  c  >  1  N  satisfies 

| Win,  ml |  '  A  a  C ^ 

| Win,  ml  -  Win  -  1,  m  -  1)|  *  A*  a 

provided  nA  -  t  satisfies  t  <  S,  whete  A  -  Mt.  C,  ,  C,,  v-  ,  v"  ,  i',  ,  A)  '  0 

1  I  J  4  9 

bounded  away  from  sero,  for  fixed  a  '  1 ,  C, . l'  Independently  v'f  A  * 

1  n 

Proof •  You  first  establish  the  representation 
n-1  _  n-1  _ 

win,  ml  •  \  W(l,  ntm-  2)  -  1)  -  \  W|t),  nrm-  21)  ♦  }'  A  Q(k,  1) 

1-0  1-1  Cln.wl 

n  '  1 

where  C(n,  ml  -  l(K,  1)  i  0  *  k  »  n ,  m  -  n  «  k  «  1  '  m  *  n  V)  and  yMO,  mt  * 
Oil,  m)  0. 


The  representation  (6.6)  clearly  holds  for  n-1.  Assume  that  it  holds  for 


W(k,  m) ,  m  f  S  and  k  <  n.  Write  the  difference  equation  in  the  form 

5(n,  m)  »  W(n  -  1,  m  +  1)  ♦  W(n  -  1,  m  -  1)  -  W(n  -  2,  m) 

2  -  2  -  - 
+  A  V  (n  -  l)W(n  -  1,  m)  +  A  W(n  -  1,  m)D  W(n  -  1,  n  -  1) 

2  -  .  -  3 

-  A  W(n  -  1,  m)D  W(n  -  1,  n  -  1)  -  A  P(n,  m) 

Note  that  (C(n  -  1,  m  +  1)  u  C(n  -  1,  m  -  11)  \  C(n  -  2,  m)  -  C(n,  m)  \  {(n,  m)}. 

Also, 

n-2  _  n-1 

l  5(1,  n  +  m  -  2J  -  1)  -  l  5(0,  n  +  m  -  2 j ) 
j-0  j-1 

n-2  n-1 

+  l  5(1,  n  +  m  -  2j  -  3)  -  £  5(0,  n  +  m  -  2j  -  2) 

j-0  j-1 

(6.7) 

n-3  n-3 

-  I  5(1,  n  +  m  -  2j  -  3)  +  ]*  5(0,  n  *  m  -  2j  -  2) 

j-0  j-1 

n  n 

-  £  5(1,  n  m  -  2j  -  1)  -  £  W  (0 ,  n  +  m  -  2j) 

j-0  j-1 

Obtain  from  (6.4),  (6.7)  and  the  induction  hypothesis 

n  n 

5(n,  m)  -  ][  5(1,  n  +  m  -  2j  -  1)  -  £  5(0,  n  +  m  -  2j) 
j-0  j-4 

+  l  A3  Q(k,  j)  +  A2  V (n  -  l)5(n  -  1,  m) 

C(n,m)\{  (n,m) } 

+  A2  W(n  -  1,  m)D  5(n  -  1,  n  -  1)  -  A*  5(n  -  1,  m)b  5(n  -  1,  n  -  1)  -  A3  T(n,  m) 
The  representation  formula  therefore  holds  for  (n,  m)  if  you  set 

Q(n,  m)  =  -P(n,  m)  +  A  3{V(n  -  l)W(n  -  1,  m)  + 

(6 . 8) 

W(n  -  1,  m) D  5(n  -  1,  n  -  1)  -  5(n  -  1,  ra)D  5(n  -  1,  n  -  l)j 
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Denote  by  C(n)  =  sup  |o(k,  m) | .  The  next  step  is  to  estimate  Q(n)  in  terms 
k<n 
miX 

of  Q(n  -  1)  ,  C  .  u  -  1;  .  .  .  ,  5- 

-  -  -  +  - 

The  main  ingredient  is  the  estimation  of  D  W.  Recall  that  D  •  D  +  D,  ♦  Dj . 

d"  W(k,  k)  -  A_1(W(k,  k)  -  W(k  -  1,  k  -  1) )  -  A-1(W(1»  2k  -  1) 
n 

k  - 

-  W(0,  2k  -  2)  +  l  A <Q(k  -  j,  k  +  j)  +  Q(k  -  j  -  1,  k  +  j))l 

J-o 

So  )d“  W(k,  k)|  £  A  +  2k  A2  Q(k)  (6.9a) 

Similarly 

|  (D*  +  D~)W(k,  k)  |  <  A  Cx  +  2k  A2  Q(k)  (6.9b) 

These  estimates,  when  combined  with  (6.8)  and  the  hypotheses,  yield  the  required 
estimate  of  Q(n) . 

The  remainder  of  the  proof  is  a  finite  induction  argument. 

Claims:  Suppose 

f  -1  C1 

x  =  nA  <_  min  (  (8(C3  +  A  k  CQ) )  ,  (*  -  1) 

Then 

(i)  |Q  (n,  m)|  <  C  5  2(C5  t  2  Cj  C3  I  K  (CQ  C4  +  2  A  CQ  C^J 

(ii)  |w(n,  m)  |  <1:C0 

(iii)  |w(n,  m)  -  W(n  -  1,  m  -  1)  |  <_  A  k  Cj 

For  n  =  1,  (i)  holds  since  (2(1,  •)  «  0  t  (ii)  and  (iii)  are  just  the  hypotheses  of 
the  Lemma.  Suppose  (i)  ,  (ii),  (iii)  hold  for  n  <_  n  -  1.  Then  from  (6.8),  (6.9)  you 
obtain 

|(2(S,  m)|^Cj  +A_1(C4  A  k  Cq  +  2  C3  A  C3  +  4  C3  A2(n-  1)0  + A  k  Cq(2  A  +  4(S  -  1)A‘  0)  ) 

<C  +  K  C  C  *  2  C,  C,  H  C,  A  ii  j  H  K  C.(2  C,  H  n  A  Q] 

—  5  04  13  3  Ul 
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I 


-  C_  +  2  C.  C  +  k  (Cn  C  +  2  A  C  C  )  +  4  n  A(C  +  A  .c  C  )Q 
5  1304  01  3  0 

The  hypothesis  of  the  claim  is  that  x  =  nA  ^  (8(C^  +  A  k  Cq)  )  *  among  other  things. 
Together  with  the  definition  (i)  of  Q,  this  implies  that  the  above  is 

ijQ+fc=Q 

which  establishes  (i) . 

The  second  part  of  the  hypothesis  on  x,  namely 


means  exactly  that 

0  -  2  +  x  Ci  '  (|C  '  1)C0  (6.10) 

On  the  other  hand,  the  representation  (6.6) ,  estimated  according  to  (i)  and  the 
hypotheses  of  the  Lemma,  gives 

n-1 

| W(n,  m) |  =  |  l  {W(l,  n  +  m  -  2j  -  1)  -  W(0,  n  +  m  -  2j  -  2)} 

3=0 

-2  3 

+  W(0,  m  -  n)  +  l  a3  Q(k,  j)|  <  S  A2  ^  +  4  C  +  =  A(Cq  +  X  C  +  X2  £) 

C(n,m) 

which  is 

iiKC0 

if  and  only  if  (6.10)  holds.  Finally,  using  (6.6)  as  before,  you  show  that 

|w(n,  m)  -  W(n  -  1,  m  -  1)  |  <  A2  C1  +  2  A3  n  Q  =  A2  (1  +  2  A  ii  £-) 

Q  - 

<  A2  C,(l  t  2(k  -  1)  4-  f)  =  A2  sC.  . 

1  2Q.  1 

The  claim  is  thus  proven,  and  with  it  the  Lemma,  if  you  take 
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Set 


Qo  -  2  c5  +  2  ci  S  +  CQ(C4  ♦  2  Cl) 


6.  *  M*'  Cn'  Ci'  Co'  C„'  CJ 
0  0  012345 


,  f  i  *  -  i)ci  r/fiV 

’lnl8(C3  +  C0)  '  Q0  '  IAqJ 


2  (k  -  1)C, 


1/2 


and  define  o  =  T/6  .  Suppose  that  A  is  so  small  that  A  k  <  1.  Then  for 


1  £  j  <  o  certainly 


Also 


8(C3  +  A  <J  CQ)  <  8(C3  +  CQ) 


Qj  =  Q(K  ,  ic  ^  CQ.  K-’  C1#  C2>  C3,  C4,  C5)  =  2  C5  +  2  k  j  C3  C3  +  (<  j  CQ  C4  +  2  A  K  2’  CQ  Cj) 


<2C  +2KJ  C.  C,  +  K  •)+1  C„  (C  +  2  C,  )  <  K3  (2  Cc  ♦  2  C,  C,  +  <  C  (C  +  2  C,  D  =  Qn 

~5  13  04  1—  5  13  04  1  0 

Hence,  if 

v  2  ,  ,.„-,l/2 


that  is 


then 


,\(h\  ..xiar 

_l\2o  /  20  J  20 


2 

0<T  yUCj  -  («-  1)CQ 


0  <  T2  K-1  +  T  ^  Cj  •  (K  •  1)  1  CQ 

2  j  j 

<  T  -J-+TICJC1-OC-l);lC0 


hence 


'f-ifiV  .  2  — ' 11  1 
V  5)  /  5l 


,ic 
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Referring  to  the  definition  (6.11)  of  6,  you  see  that,  6^  >  6^,  j  =  1,  .  .  .  ,  o. 
So  you  can  replace  6^  by  6^  in  the  estimates  at  the  beginning  of  the  proof  to  get 

[  W  (.n ,  m)  |  <_  A  K  3  CQ 

for  (j  -  Dig  nA  <_  j  6  ,  j  =  1,  .  .  .  ,  o.  For  j  =  o,  this  gives  the  estimate 

we  want  for  W.  The  second  estimate  now  follows  in  exactly  the  same  way. 

q.e.d. 

Theorem  6.3.  Suppose  W  satisfies  (6.4),  (6.5).  Then  for  (n,  m)  e  C*,  arbitrary 


<  >  1,  and  A  sufficiently  small. 


|W(n,  m)  <  A  k  C„ 


|6  W(n,  »)  |  <  4  k  C 

where  C  ,  are  entire  functions  of  «  \  A,  T,  II  fil  ,  II D'  Fll  ,  II D+  D  Fll  and  o  >  0 
is  the  minimum  envelope  of  entire  functions  of  these  arguments. 

A 

Proof.  First  observe  that  the  difference  scheme  (6.4)  restricts  to  C  .  Also,  all 
of  the  constructions  used  in  the  proofs  of  Lemma  6.1  and  Corollary  6.2  restrict 

T 

Now  (6.5)  can  be  written 

»  (1,  m)  *  A^  B(m,  A) 

It  follows  that  the  hypotheses  of  Lemma  6.1  are  satisfied,  with  estimated  by  the 

bound  for  B  mentioned  at  the  end  of  section  5,  =  A  C^,  C,.  estimated  by  the 

bound  for  P  given  in  section  5,  and  and  estimated  by  the  bounds  for  W, 


V  given  in  section  4.  Note  that 


D  W(n,  n)  =  -(D  +  D  +  D  ) W (n,  n) 

n  2  1 


-A  {W(n,  n)  -  W(n  -  1,  n  -  1)  +  W(n,  n  +  1)  -  W(n,  n) 

+  W(n,  n)  -  W(n  -  1,  n) } 


-A  (W(n,  n)  -  W(n  -  1,  n  -  1)  +  W(n,  n  +  1)  -  W(n  -  1,  n)  } 
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and  both  of  the  differences  in  this  expression  are  estimated  by  Corollary  b.2  which 


leads  to  the  second  inequality  in  the  statement.  Note  that  the  definition  of  5 


in  the  proof  of  Corollary  6.2  shows  that  o  is  as  described 


Corollary  6.4.  For  the  solution  W  ,  V  of  (6.1)  and  the  solution  W,  V  of  (5.14) 


we  have 


We  now  turn  to  the  second-order  consistent  scheme  mentioned  in  section  2,  which 


we  reproduce  here 


As  it  stands,  this  scheme  is  incomplete:  the  value  V(0)  is  not  determined 


from  the  data.  In  order  to  preserve  second-order  consistency  we  add  the  following 


one-sided  difference  approximation  to  V(0) ,  obtained  from  the  GL  equation  and  the 


boundary  condition  (2.2c) 


Theorem  6.5.  For  any  p  >  0,  there  exists  d  >  0  so  that  for  0  <  d  <  A  ,  the 


solution  W 


Proof.  The  idea  is  to  compare  W  and  W  in  the  same  way  W  and  W  were 


compared  above 


Write 


Then  (6.1)  can  be  rewritten 


V (C) A  )W  (0,  m) 


V (0) A  -  HA)W  (0,  m)  +  (1 


P  (n,  m.  A)  =  A"  (D*  W  (n,  n)  -  D  W  (n,  n))W  (n,  m) 


A  (D,  W  (n,  n  +  1)  -  D  W  (n,  n))W  (n,  m) 


D~  W  (n,  n  +  1)  +  D  W  (n,  n  +  1)  -  D  W  (n,  n))W  (n,  m) 


D,  W  (n,  n))W  <n,  m) 


According  to  Proposition  4.3(c),  P  can  be  estimated  by  an  entire  function  of 


T,  A,  and  norms  of  F  and  its  first  and  second  differences 


W,  we  obtain 


_  _ 


From  this  point,  the  proof  follows  exactly  the  proofs  of  Lemma  6.1,  Corollary  6 


Theorem  6.3,  and  Corollary  6.4,  and  we  leave  the  details  to  the  reader.  We  note  in 


analogous  to  those  of  Theorem  6.3 


passing  that  estimates  appear  for  W'  and  D'  W 


Remark.  It  seems  likely  that  similar  results  Could  be  obtained  for  higher-order 


consistent  schemes  for  the  Chudov  problem 


Results  of  Numerical  Experiments 


wv  present  the  results  of  two  series  of  numerical  experiments  carried  out  at  the 


University  of  Wisconsin  -  Madison  MACC  facility.  The  computations  were  performed  in 


single  precision  on  the  UNIVAC  1110 


Both  series  of  experiments  involved  a  FORTRAN  program,  detailed  below,  which 


implements  the  second-order  Chudov  Scheme  (6.12),  In  both  cases  V  is  to  be 


computed  on  [0,  11,  i.e.  T  »  1.  The  series  differ  in  the  mode  of  data  generation 


In  the  first  series,  the  inverse  problem  for  F(t)  F  (const.)  is  solved  by  way 


of  (6.1 2).  The  approximate  potential  V(n)  =  V(nA)  is  compared  with  the  exact 


potential  VE(n)  »  VE(nA),  which  is  known  in  closed  analytic  form 


Both  the  maximum  (sup  norm)  error 


The  last  three  experiments  in  series  1  are  meant  to  illustrate  the  dependence 


of  the  error  on  the  lower  bound  e ,  which  for  these  examples  is  equal  to 


1,  the  exact  potential  given  by  (7.1)  clearly  admits  no 


uniform  Lipschitr  bound  in  terms  of  norms  of  F  alone  (in  fact,  the  analogous 


expression  for  W  shows  that  the  bound  (2.Sa)  is  sharp).  This  lack  of  Lipschitz 


uniformity  is  reflected  in  the  behaviour  of  the  computation .  as  predicted  by  the 


theory.  For  II  =  -.99,  the  computer  produces  garbage  for  N  *  11.  Even  for 


N  -  101,  the  errors  are  relatively  large,  although  almost  all  of  the  error  is 


concentrated  at  the  "deepest”  end  of  the  interval  (near  x  =  1 


predicted  by  the  theory) 


The  second  series  of  experiments  is  based  on  numerically  computed  F 


various  V,  H.  That  is,  VE  and  H  are  selected,  the  Chudov  problem  is  solved 


numerically  for  F  (a  version  of  the  forward  scattering  problem),  then  th,  ;  rc  ra: 


given  below  is  executed  with  various  samplings  of  this  numerically-generated  F  as 


input,  and  the  resulting  approximate  V  is  compared  with  VE.  Again,  maximum  and 


average  errors  are  displayed 


The  main  difficulty  in  this  second  series  was  the  data  (F  -)  generating 


program,  which  was  asked  in  effect  to  solve  a  discontinuous  initial  value  problem  for 


a  wave  equation.  We  finally  settled  on  a  staggered-grid  method  which  computes  F  at 


2000  points  on  the  t-axis,  which  are  then  sampled  to  produce  the  input  for  the 


inverse  program.  Despite  the  relatively  fine  grid  of  the  forward  computation,  error 


in  F  has  an  observable  effect  on  the  convergence  of  the  inversion  scheme  in 


several  examples,  in  the  (N  =  SI)  to  (N  =  101)  step 


It  is  difficult  to  assess  the  error  figures  given  in  the  tables  below,  since  the 


problem  solved  here  is  merely  a  model  problem,  with  no  physical  interpretation  what 


soever  attached.  For  the  same  reason,  we  have  not  bothered  to  predict  the  errors 


based  on  the  theory.  Both  of  these  items  would  be  interesting  for  some  of  the 


problems  mentioned  at  the  end  of  the  next  section 


<6. 12) 


1.  Read  VK (I >  ,  1-1  ,  N 

F  (1 )  ,  1-1  .  2N  -  1 

2.  H  -  F (1) 

3.  Vll)  -  2  H*  -  2  A_l(2  F(2)  -  ~  F{3)  -  ~  F(l)) 

4.  For  1*1  ,  2N  -  1 

W(l,  I)  -  F  (I ) 

5.  For  1-2  ,  2N  -  2 

W(2,  I)  -  (7  &  Vll)  -  H)WU,  1)  +  J  (Wll.  1  +  1)  +  Wtl.  I  -  1)) 

6.  For  X  -  3  >  N 

6.1  W(I,  I)  »  W(I  -  1,  I  +  1)  +  W(I  -  1,  1  -  1)  -  W(I  -  2,  I) 

+  A[W(I  -2,1-  2)W  (I  -  1,  1)1  11  +  A  W(I  -  1,  I))'1 

6.2  Vll  -  1)  -  -&'1(W(I,  I)  -  W(I  -  2,  I  -  2)) 

6.3  For  J-I+l  ,  M-I+l 

6.3.1  W(l,  J)  -  WU  -  1,  J  +  1)  +  W(I  -  1,  3  -  1)  -  W(I  -  2,  J) 

+  A*  V{I  -  1)W(I  -  1,  J) 

?.  MAX  ERROR  »  MAX  |  V  (I )  -  VE(I)| 

1<I,  <N-1 


8.  AVE  ERROR  -  IN  -  l)'1  £  | V (I )  -  VE(I) 

1-1 
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*  'Ait  / 


n  mW?!'?  •••*"*  ’  |K7^7l>KF»«*^,iPWWjHW» '(■•■  •  ft1' .-  •  ••>«;■■-»••  ~yiPgM»«a*isjgi  |  H  i  ;nym  <f*!)l>ll".  |  'HWftiWf! 


A. 


B. 


c. 


D. 


E. 


tabu;  ii 


V  (X)  0  H  -  2 


N 

A 

MAX. 

ERROR 

AVG. 

ERROR 

11 

.1 

.1 

.97 

-1 

21 

.05 

.30 

-1 

.18 

-1 

51 

.02 

.53 

-2 

.15 

-2 

101 

.01 

.18 

-2 

.57 

-3 

V(X) 

£  0 

H  -  -2 

11 

.1 

.9 

-1 

.35 

-1 

21 

.05 

.24 

-1 

.11 

-1 

51 

.02 

.33 

-2 

.29 

-2 

101 

.01 

.15 

-2 

.14 

-2 

V(X)  -  1 

-  X 

M  -  0 

11 

.1 

.25 

-2 

.17 

-2 

21 

.05 

.61 

-3 

.41 

-3 

51 

.02 

.10 

-3 

.6 

-4 

101 

.01 

.30 

-4 

.1 

-4 

V(X) 

-  X(1  -  X) 

H  «  0 

11 

.1 

.42 

-2 

.36 

-2 

21 

.05 

.11 

-2 

.95 

-3 

51 

.02 

.26 

-3 

.19 

-3 

101 

.01 

.15 

-3 

.90 

-4 

V(X) 

-  COS  (2X) 

II  -  0 

11 

.1 

.64 

-2 

.37 

-2 

21 

.05 

.17 

-2 

.93 

-3 

51 

.02 

.27 

-3 

.15 

-3 

101 

.01 

.7 

-4 

.5 

-4 

F 


V(X)  -  5  COS(IOX) 


II  -  0 


58.  Conclusion 


We  have  still  to  discuss  the  relation  of  our  results  to  work  of  other  author-, 
and  the  relative  efficiency  of  our  method  for  solving  the  inverse  prol  lem.  sino 
the  first  point  bears  on  the  second,  we  begin  with  it. 

Our  results  belong  to  the  line  of  work  begun  by  Gel'fand  and  Levitan  in  111.  W. 
have  discussed  the  relation  of  previous  work  in  this  line  to  our  results  for  the 
continuum  inverse  problem  in  111,  and  we  defer  further  discussion  of  the  various 
approaches  to  the  continuum  problem  to  another  place.  We  restrict  our  discussion 
here  to  approximate  methods. 

All  but  one  or  two  of  the  authors  writing  on  inverse  problems  in  the  spirit  of 

Gel'fand  and  Levitan  base  their  work  on  certan  a  linear  integral  equation.  This  is 

either  the  Gel’ fand-Levitan  integral  equation  or  the  Marchenko  integral  equation, 
depending  on  whether  the  incident  waves  are  incident  at  the  origin  or  at  infinity 
(the  latter  is  always  the  case  for  tile  whole-line  inverse  scattering  problem)  .  Also, 
the  inverse  problems  are  usually  formulated  in  the  frequency  domain:  that  is,  the 
scattering  datum  is  either  the  spectral  function  (which  is  the  Fourier  transform  of 
our  F(t):  see  (11,  section  3),  the  (frequency-dependent)  phase  shift,  or  the 
(frequency-dependent)  scattering  matrix  or  reflection  coefficient. 

Among  the  authors  whose  work  on  approximate  solutions  to  inverse  problems  fits 
into  the  framework  just  described  are  Case  (9),  Case  and  Kac  110],  Case  and  Chiu  111), 

Ware  and  Aki  (121,  Berryman  |13),  and  Green  and  Berryman  (141. 

P.  Gopillaud  |151  has  given  a  slightly  different  (and  faster)  algorithm  for  the 
normal  incidence  elastic  waves  inverse  problem  for  layered  media.  Ware  and  Aki  1 1 1  , 
Greene  and  Berryman  (13),  and  Berryman  (141  develop  this  idea  further  and  show  that 
the  Gopillaud  and  discrete  Gel ' fand-Levitan-Marchenko  approaches  are  equivalent,  in 
various  senses. 


-40- 


Our  work  differs  from  all  of  the  above  mentioned  work  in  many  respects,  the 
following  three  of  which  we  believe  most  important: 

1)  We  base  all  of  our  results  on  the  nonlinear  Volterra  equation  (GL  in  section  2) 
rather  than  the  linear  Gel 1  f  and- Levi  tan  or  Marchenko  equations.  Though  these  are,  in 
principle,  equivalent,  it  seems  difficult  to  extract  the  proper  stability  results  from 
the  linear  equations. 

2)  Both  the  formulation  of  our  inverse  problem  and  its  solution  are  time- 
dependent.  It  is  clear  from  the  results  of  [1)  that  at  least  some  of  the  "frightful 
instability"  attributed  by  Wheeler  (161  to  the  Gel ' fand-Levitan  method  enters  by  way 
of  the  conditionally  convergent  Fourier  integrals  required  to  pass  between  the  time- 
domain  and  frequency-domain  problems. 

3)  Most  important  for  application  and  generalization  is  numerical  stability. 

Some  stability  result  is  required  to  guarantee  that  each  algorithm  will  actually 
converge  as  &  ■*  0  to  the  solution  of  the  continuum  problems.  None  of  the  above 
authors  seem  to  provide  the  necessary  estimates  to  conclude  convergence.  To  provide 
such  estimates  is  exactly  the  point  of  the  present  work. 

The  Cholesky  decomposition  estimates  of  section  4  should  also  provide  estimates 
for  the  solution  of  the  discrete  versions  of  the  linear  Gel ' fand-Levitan  integral 
equation,  hence  imply  convergence  of  the  various  algorithms  so  constructed.  The  same 
should  be  true,  with  proper  attention  to  the  behaviour  of  the  potential  at  infinity, 
for  the  discrete  Marchenko  equation  approaches. 

We  conjecture  that  the  Gopillaud  algorithm  is  actually  closely  related  to  our 
Chudov  scheme,  hence  should  inherit  the  stability  properties  derived  here. 

The  approaches  to  the  inverse  problem  based  on  the  various  integral  equations 

4  3 

have  costs  proportional  to  N  or  N' ,  depending  on  implementation  (N  =  number  of 

linear  gridpoints,  as  usual).  The  Gopillaud  scheme  (and  modifications  -  1121,  1131, 

2 

and  [141)  has  cost  proportional  to  N  .  The  Chudov  schemes  investigated  here  also 

2 

have  cost  proportional  to  N  ,  hence  seem  to  be  optimally  cheap  amongst  "exact" 


inversion  methods. 


”  — - 


In  applications  of  inverse  problems,  for  instance  in  exploration  geophysics, 
physical  chemistry,  and  ultrasound  tomography,  various  ''approximate"  inversion  methods 
are  used.  Indeed  these  problems  involve  more  than  one  space  dimension,  and  "exact" 

Gel 1 fand-Levitan-Marchenko  methods  have  yet  to  be  extended  in  any  useful  way  to  higher¬ 
dimensional  problems.  These  approximate  methods  are  based  on  well-known  approximate 
solution  methods  for  the  relevant  partial  differential  equations,  for  instance  Born 
series  approximation  (formal  perturbation  series)  or  JWKB  methods  (geometric  optics). 
The  approximate  solutions  are  then  inverted  to  obtain  (one  hopes)  approximate  solu¬ 
tions  to  inverse  scattering  problems.  One  can  also  apply  these  methods,  as  detailed 

for  instance  in  (17),  Chapter  XVI,  to  the  simple  inverse  problem  of  this  paper.  One 

2 

obtains  schemes  with  cost  proportional  to  N  (or,  for  the  Born  approximation, 

N  log  N  if  one  employs  the  Fast  Fourier  Transform) .  Our  Chudov  algorithm  therefore 
has  more  or  less  the  same  cost  as  these  "approximate"  methods.  Of  course,  the  Chudov 
algorithms  also  enjoy  stability  properties,  as  we  have  shown,  which  imply  the 
possibility  of  rigorous  error  estimation.  No  such  possibility  is  available  for  the 
"approximate"  methods,  to  our  knowledge. 

We  describe  in  conclusion  some  other  problems  to  which  our  methods  apply.  The 
analytical  details  (results  analogous  to  Theorems  1  and  2,  and  estimates  (2.5a,  b,  c)) 
have  been  worked  out  in  (18)  for  an  inverse  problem  for  the  acoustic  wave  operator 


a2  2  a2 

°C  “  2  *  C  x>  2 

at  ax 

and  in  [19)  for  some  inverse  problems  for  general  hyperbolic  systems  in  two  variables 
of  first  order,  with  constant  sound  speeds.  Numerical  results  should  also  follow  by 
the  techniques  of  this  paper.  In  these  problems,  the  incident  waves  are  incident  at 
some  finite  point  in  the  (one-dimensional)  medium.  Inverse  scattering  problems,  in 
which  waves  are  incident  at  infinity,  should  also  submit  to  roughly  the  same  approach, 
although  the  technical  details  remain  to  be  worked  out.  We  point  out  that  inverse 
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problems  for  the  wave  operator  1  typically  are  approached  through  dependent-  and 
independent- variable  transformations  which  convert  it  into  the  potential  perturbation 
of  0  treated  in  this  paper  (112),  1141) .  Neither  this  possibility  nor  the 
Gel ' fand-Levitan-Marchenko  linear  integral  equation  approach  are  available  for  the 


general  hyperbolic  first-order  system  with  more  than  one  variable  (unknown)  sound 
speed.  We  have  succeeded  in  solving  some  inverse  problems  for  such  systems  by  combin 


ing  the  techniques  of  1 18]  and  119) 


Finally  we  mention  that  the  artifacts  of  our  approach  to  inverse  problems 
(GL  equation,  Chudov  boundary- value  problem)  are  also  available  for  higher-dimensional 


inverse  problems.  We  have  not  yet  derived  the  necessary  a  priori  estimates  to  proceed 


with  such  an  extension  of  the  theory,  however 
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